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Abstract 

We perform large-scale simulations of directed sandpile models with both 
deterministic and stochastic toppling rules. Our results show the existence 
of two distinct universality classes. We also provide numerical simulations of 
directed models in the presence of bulk dissipation. The numerical results 
indicate that the way in which dissipation is implemented is irrelevant for 
the determination of the critical behavior. The analysis of the self-affine 
properties of avalanches shows the existence of a subset of super-universal 
exponents, whose value is independent of the universality class. This feature 
is accounted for by means of a phenomenological description of the energy 
balance condition in these models. 

PACS numbers: 05.65. +b, 05.70.Ln 
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I. INTRODUCTION 



Sandpile cellular automata are the most famous example of self-organized critical (SOC) 
behavior . Under an external drive consisting of a slow addition of sand (energy) grains 
and the action of dissipation through the loss of energy on the lattice boundaries, these 
models reach a stationary steady state. In the limit of infinitesimal driving and dissipation 
(this last achieved in the thermodynamic limit), the stationary state of sandpile models 
exhibits diverging response functions associated to a characteristic avalanche dynamics. This 
is the hallmark of a critical behavior that has attracted an enormous amount of interest as 
a plausible explanation of the avalanche-like critical behavior empirically observed in many 
natural systems . 

Sandpile models have been at the center of an intense research activity made of both an- 
alytical studies and numerical simulations. Despite the simple definition of these automata, 
it turns out that their full analytical understanding is a very problematic task |Q. As a fur- 
ther complication, also the numerical inspection of these models results to be particularly 
difficult. For example, the precise identification of universality classes has resisted for many 
years even the most careful numerical analysis, and only recent results have partially settled 
this problem On the other hand, these refined analyses have pointed out that several 

sandpile models do not follow the simple finite size scaling (FSS) form usually adopted in 
the description of critical behavior [Q. For instance, the more sophisticated multiscaling 
approach [|TU]-[T^ seems to be required for a full description of the scaling properties of the 
original Bak, Tang, and Wiesenfeld (BTW) model IHQ]. 

Many sandpile features have been underlined as the possible origin of these scaling 
anomalies. The deterministic dynamical rules of the BTW model induce nonergodic effects 
P], that are certainly missing in stochastic models, such as the Manna model 1131,1^, which 
shows a perfect FSS behavior, even for moderate system sizes. A further complication of 
sandpile automata stems from the peculiar role of the boundary dissipation, that makes the 
lattice size scaling entangled with the system dynamics. In such cases, the thermodynamic 
limit is essential for the dissipative dynamics of large avalanches. A clear understanding of 
the interplay between dissipation and size scaling has not yet been achieved and it has been 



recently the subject of several studies |11,14 



In this paper we address some of the aforementioned problems in the case of directed 
sandpile models [p!0| , p!5| -[T8[| . In this case Dhar and Ramaswamy obtained an exact solution 
for the Abelian deterministic directed sandpile (DDS) [|l^], that can be used as a milestone 
to check the numerical simulation analysis. Directed sandpiles thus become an interesting 
test field to study how the critical behavior is affected by the introduction of stochastic 
elements and dissipation. We perform large scale numerical simulations of two directed 
sandpile automata: the deterministic directed sandpile model |15| and the stochastic directed 
sandpile model |]TB[. We study both models in the case of boundary and bulk dissipation 
||19|^2[. We find, in agreement with the results in Ref. |]18[, that the models define two 



different universality classes. In addition we show that the universality class of the models 
does not depend on the way in which dissipation is implemented. Finally we analyze the 
properties of anisotropic models in which the dynamics is not fully directed In 
this case we observe that on large scales the critical behavior is the same of that of fully 
directed models. Results for the stochastic model are compared with a recent theoretical 
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approach by Paczuski and Bassler p4|, that provides values for the critical exponents in 



perfect agreement with numerical simulations. These results are also recovered in Ref. ^5 



The numerical analysis also points out that some critical exponent values, such as the 
correlation length exponents or the affinity exponent (to be defined later on), are independent 
of the particular universality classes and common to all models considered. In order to 
explain this numerical evidence, we provide a phenomenological characterization of directed 
sandpiles based on the basic symmetries introduced by the conserved dynamics of these 
automata. Following balance of energy arguments inspired in Refs. [p^ -p8[|, we derive a 
series of results and predictions on the value of critical exponents which are a straightforward 
consequence of conservation. These general results can be considered as super-universal, 
because they characterize the critical behavior of all directed sandpiles with local dynamical 
rules, independently on the specific universality class. The results presented here provide 
a general picture of directed models and the role of boundary and bulk dissipation in the 
process of self-organization. 

The paper is arranged as follows. In Sec. II we introduce and define the various directed 
models considered. Sees. Ill and IV present and discuss, from the standpoint of universality, 
the numerical results for directed models with boundary and bulk dissipation. In Sec. V 
we introduce anisotropic models, and present the numerical results obtained, in comparison 
with those of directed models. Sec. VI is devoted to an analytical approach based on the 
conservation of energy. Finally, in Sec. VII we draw our conclusions and perspectives. 



II. DIRECTED MODELS 

Sandpile models are usually defined on a ci-dimensional hyper-cubic lattice of size L. To 
each node of the lattice is assigned an integer variable Zi, called "energy". Energy is added 
to the system uniformly at randomly chosen sites {zi Zi + 1). When a site becomes active, 
that is, when its energy becomes larger than or equal to a certain threshold Zc, it topples. 
A toppling site loses an energy Zc, that is distributed among its neighbors according to a 
certain set of rules. The neighbors that receive energy can become active and topple on 
their turn, thus generating an avalanche. The slow driving condition is effectively imposed 
by stopping the random energy addition during the avalanche spreading. This means that 
the driving time scale is infinitely large with respect to the toppling characteristic time scale. 

The models we consider in this Section are directed, in the sense that the energy is 
always transported along a preferred fixed direction. We denote this preferred direction 
by the coordinate x\\, whose positive direction is usually defined as "downwards". The 
transverse direction (subspace of dimension d — 1 perpendicular to x\\ ) will be denoted by 
x±. 

The toppling rules of the models define two main classes: 

(i) Deterministic directed sandpile (DDS): In d dimensions, the threshold is set to Zc = 
2d + 1. When a site in a given hyper-plane xy topples, it sends deterministically one grain 
of energy to each one of its nearest and next-nearest neighbors on the hyper-plane x\\ + 1 
(see Fig. |l]a)). Our definition is somewhat different from the original model of Dhar and 
Ramaswamy |TB|, in both the driving and the orientation of the lattice. Both models, 
however, are expected to share the same universality class, being deterministic and directed. 
Numerical simulations confirm indeed this point []T^|. 
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FIG. 1. Toppling rules in d = 2 for directed sandpiles. Filled circles represent active (toppling) 
sites; empty circles are stable sites. In the deterministic model (a) an active site sends one grain 
to each of its three neighbors on the next downwards row. In the stochastic models, exclusive (b) 
and nonexclusive (c), one grain is sent to two randomly chosen downwards neighbors. 



(ii) Stochastic directed sandpile (SDS): In this case, the threshold is Zc = 2, independently 
of the dimensionality of the lattice. When a site in the hyper-plane xy topples, it sends two 
grains of energy to two sites, randomly chosen among its 2d — 1 nearest and next-nearest 
neighbors on the hyper-plane -|- 1. The toppling rules of this model can be defined 
exclusive if the two energy grains are always distributed on different sites. Fig. |l]b). On 
the other hand, the model can be defined non- exclusive if the dynamics allows the transfer 
of two energy grains onto the same site. Fig. |l|c). We therefore report simulations on the 
exclusive stochastic directed sandpile (ESDS) and on the non-exclusive stochastic directed 
sandpile (NESDS). In spite of the stochastic nature of these models, we must bear in mind 
that they are nevertheless Abelian [Q. The discussion therefore focuses on the difference 
between stochastic and deterministic models. 

Once the toppling rules have been determined, the models are finally defined by specifying 
the dissipation mechanism. For systems with boundary dissipation, we impose periodic 
boundary conditions in the transverse directions x± and open at the hyper-plane x\\ = L. In 
this way, the models are locally conserved; energy can only leave the system at the bottom 
of the lattice. In models with bulk dissipation, we impose periodic boundary conditions in 
both the X|| and x± directions. Dissipation is implemented by allowing a toppling site to 
loose an energy without transferring it with probability p [^,^. This means that, on 
average, an energy e = z^p is dissipated in each toppling. In the limit e — 0, the system 
shows critical behavior . 



In the stationary state we can define the probability that the addition of a single energy 
grain is followed by an avalanche of toppling events. Avalanches are then characterized by 
the total number of topplings s and the time duration t. In the limit of infinitesimal driving 
(slow driving condition) the system shows scaling behavior and the probability distributions 
of these quantities follow the finite-size scaling (FSS) forms: 



Pis) 
Pit) 



s-^'Gis/s, 



(1) 

(2) 



where Sc and tc are the characteristic size and time, respectively. The exponents and 
Tt characterize the critical behavior and define the universality classes to which the models 
belong. In the critical region the characteristic time and size are determined only by the 
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system size L or the dissipation e, in the case of boundary and bulk dissipation, respectively. 
In directed models, the affinity exponent ( is of particular importance; it relates the avalanche 
characteristic lengths in the perpendicular direction, C,±, and in the parallel direction, ^y, 
through the relation ^± ~ ^jj. This exponent characterizes the degree of anisotropy due to 
the preferential direction present in the transport of the energy. In other words, it expresses 
the self-affine properties in the scaling of avalanches. A general result concerns the average 
avalanche size (s), that also scales linearly with L [Il0| , p!5| , p3| : a new injected grain of energy 
has to travel, on average, a distance of order L before reaching the boundary. In the 
stationary state, to each energy grain drop must correspond, on average, an energy grain 
flowing out of the system. This implies that the average avalanche size corresponds to the 
average number of topplings needed for a grain to reach the boundary; i.e., (s) ~ L. The 
same result can be exactly obtained by inspecting the conservation symmetry of the model 
as we shall see in Sec. |VI[ 

For the DDS, the exact analytical solution in = 2 yields the exponents = 4/3 and 
Tf = D = 3/2 |T^. The upper critical dimension is found to be = 3, and it is also possible 
to find exactly the logarithmic corrections to scaling [15,^|. The introduction of stochastic 
ingredients in the toppling dynamics of directed sandpiles has been studied recently in a 
model that randomly stores energy on each toppling |T^. This model is strictly related 
to directed percolation and defines a universality class "per se". In our case stochasticity 
affects only the partition of energy during topplings, and there is a priori no obvious relations 
between the critical behavior of these models. 



III. NUMERICAL SIMULATIONS WITH BOUNDARY DISSIPATION 

In this Section we report results from computer simulations of deterministic and stochas- 
tic directed sandpiles, performed with boundary dissipation. The system sizes considered 
range from L = 100 to L = 6400. The statistical distribution functions have been computed 
averaging over 10'' nonzero avalanches. 

In the case of boundary dissipation, the lattice size L is the only characteristic length 
present in the system. Approaching the thermodynamic limit {L oo), the avalanche 
characteristic size and time in Eqs. |1| and ^ diverge as Sc ~ and tc ^ L^, respectively. The 
exponent D defines the fractal dimension of the avalanche cluster and z is the usual dynamic 
critical exponent. The directed nature of the model introduces a drastic simplification, since 
it imposes z = 1. In order to compute the different exponents characterizing the dynamics of 
the avalanches, we have performed the moment analysis of the distributions, in analogy to the 
method developed by De Menech et al. [lllJl^ . We define the g-th moment of the avalanche 



size distribution on a lattice of size L as (s^)^ = / dss^ P{s). If the FSS hypothesis (P is 
valid in the asymptotic limit of large s, then the g-th moment has the following dependence 
on system size: 

(s")^ = j dyy^'i-^^^ G{y) ~ L'^^^'?). (3) 

The exponent ^^(g) = D{q + 1 — r^) is computed as the slope of the log-log plot of (s'')^ as 
a function of L. For large enough values of q (i.e., away from the region where the integral 
in is dominated by its lower cut-off), one can compute the fractal dimension D as the 
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slope of (Tsiq) as a function of g: D = das{q)/dq. On the other hand, as we have argued in 
the previous Section, the first moment must scale linearly with L, which imposes o"s(l) = 1- 
Once D is known we can estimate using the relation (7^(1) = D{2 — r^) 

Along the same lines we can obtain the moments of the avalanche time distribution. In 
this case, {t'i)L ~ L'''^'^\ with dat{q)/dq = z. Analo gous considerations for small q apply 
also for the time moment analysis. Here, an estimate of the asymptotic convergence of the 
numerical results is provided by the constraint z = 1, that must hold for large enough sizes. 
Then, the Tt exponent can be found using the scaling relation (2 — r^) = <7t{l). 

Once the exponents have been estimated numerically, we can check the accuracy of the 
moment analysis' predictions using the FSS hypothesis. If the FSS hypothesis of Eq.s (|l],|D 
is correct, then the plots of the distributions, under the rescaling s s/L^ and P{s) — >■ 
P{s)L^'^'' and correspondingly t t/L^ and P(t) — > P(t)L^'^\ should collapse onto the same 
universal function, for different values of L. 

In Table | we report the exponents found for the DDS, ESDS, and NESDS models in 
d = 2. Figure ^ shows the moments crs(g) and (Jt{q). Figures | and || plot the FSS data 
collapse for sizes and times, respectively. The exponents obtained for the DDS are in perfect 
agreement with the expected analytical results. This fact supports the idea that the system 
sizes used in the present work allow to recover the correct asymptotic behavior. Results 
for the ESDS and NESDS are identical within the error bars, pointing out that these two 
models are in the same universality class. On the other hand, the obtained exponents prove 



Model Ts D Tt z C 

DR 1/3 3/2 3/2 1 1/2 

DDS 1.34(1) 1.51(1) 1.51(1) 1.00(1) 0.50(1) 

ESDS 1.43(1) 1.74(1) 1.71(3) 0.99(1) 0.51(1) 

NESDS 1.43(1) 1.75(1) 1.74(4) 0.99(1) 0.51(1) 



TABLE I. Critical exponents for directed sandpiles with boundary dissipation in d = 2. DR: 
Dhar and Ramaswamy's exact solution; DDS, deterministic directed model; ESDS and NESDS, 
stochastic directed models. Figures in parenthesis denote statistical uncertainties. 
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FIG. 3. Data collapse analysis of the integrated avalanche size distribution for the d = 2 
stochastic models with boundary dissipation a) ESDS and b) NESDS. System sizes are 
L = 400, 800, 1600, 3200, and 6400. 
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FIG. 4. Data collapse analysis of the integrated avalanche time distribution for the d = 2 
stochastic models with boundary dissipation a) ESDS and b) NESDS. System sizes are 
L = 400, 800, 1600, 3200, and 6400. 



beyond any doubts that deterministic and stochastic directed sandpile models do not belong 
to the same universality class. 

We have also directly computed the characteristic lengths in the parallel and transversal 
directions, and as a function of the system size. The anisotropy of the system is 
reflected in the different definitions of both characteristic lengths. In this sense, we define 
them with the same spirit as in directed percolation [Q. 

Consider a given avalanche, labeled a, that has started at the site {x^^\x^^^), and has 

affected the set of different sites {(xy^^ x^±)}, for z = ... a — 1 (i.e., it has covered an area 
a). Let us define the quantities 

= -EVr -4^1 (4) 

^ 1=1 



and 
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logio L 

FIG. 5. Correlation lengths and as a function of L for the models with boundary dissipa- 
tion DDS (o), ESDS (A), and NESDS (o). The dashed lines are guides to the eye with slope 1.00 
and 0.50. 
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Furthermore, let us define R\\{a) and I^^{a) as the averages of the previous quantities, over 
all avalanches of the same fixed area a. Let be P{a) the probability of observing an avalanche 
of area a. We define the correlation lengths by 



Eg R\\{a)aP{a) 



(6) 



The different definitions and are obviously due to the different nature of the avalanche 
spreading in the directions x\\ and x±. In the former case, the spreading is isotropic, and 
thus the second moment of the relative distance distribution is needed to define a meaningful 
correlation length. In the latter case, on the other hand, the spreading is always in the 
direction of growing xy, and therefore the first moment is sufficient. 

The system being critical, both correlation lengths should scale with the system size, 
defining the exponents z/n and by 



The affinity exponent, defined by 
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(7) 



(8) 



is thus given by C = ^±/^\\- 

We have calculated the correlations lengths in the models DDS, ESDS, and NESDS, 
given by the definition (||). The results, plotted in Fig. ^ give the following dependence of 
the correlation lengths with system size for all models: 



^11 ~ L, ~ L^/^. 



(9) 
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These relations define the exponents z/y = 1 and = 1/2, and an affinity exponent ( = 1/2. 
It is interesting to note that this exponent is independent of the universality class of the 
model, defining a sort of super-universal property of directed models. 



As pointed out in Ref. [|18[, the stochastic dynamics of SDS models introduces multiple 
toppling events on the same site, which are by definition absent in the deterministic case. 
This gives rise to very a different avalanche structure, eventually reflected in the different 
asymptotic critical behavior. It is worth remarking that the universality class of SDS appears 



robust to modifications of the stochastic microscopic dynamics as pointed out in Ref. pi 
where it is shown that modifications of SDS models with stochastic toppling threshold still 
belong to the same universality class. Recently, Paczuski and Bassler ||2^, have proposed a 
theoretical approach that allows the calculation of critical exponents in directed models with 
multiple topplings. The analysis goes through the mapping of the avalanche evolution into 
the dynamics of an interface moving in a random medium, as also proposed in |l32| , |33| . This 
theoretical result gives the exponents Ts = 10/7 and Tf = D = 7 / 4, in perfect agreement 
with the values obtained by numerical simulations. Table |. The same exponent values are 



also found in the approach of Ref. [25 



IV. NUMERICAL SIMULATIONS WITH BULK DISSIPATION 

In this Section we report results from computer simulations of deterministic and stochas- 
tic sandpiles, performed with bulk dissipation. In this case, dissipation is implemented as 
described in Sec. ||. That is, in a system with periodic boundary conditions, each toppling 
site has a probability e/z^. of losing an energy Zc, and a probability 1 — e/z^ of transferring 
it to its neighbors. The dissipation rates range from e = 0.0016 to 0.0512, and the (fixed) 
system size considered is L = 6400. Statistical distribution functions have been computed 
averaging over 10'' nonzero avalanches. 

In the presence of bulk dissipation the characteristic sizes are determined by the dissi- 
pation rate e, which defines the only characteristic length in the system. Approaching the 
limit e — » 0, the avalanche characteristic size and time diverge as Sc ~ e"^" and tc ~ e""^*, 
respectively. It is also very easy to relate the mean avalanche size to the dissipation rate e. 
On average, each added grain must be dissipated in the evolution of the avalanche, resulting 
in e (s) = 1. This readily yields (s) = e~^. In this case it is extremely important that the 
characteristic length of the avalanche is always smaller than the size of the lattice used. 
This allows us to study only finite size effects introduced by the dissipation probability, 



Model Tt At ^ 

DR 4/3 3/2 3/2 1 1/2 

DDS 1.32(1) 1.50(1) 1.52(1) 1.00(1) 0.51(1) 

ESDS 1.42(1) 1.72(1) 1.70(4) 0.98(2) 0.51(1) 

NESDS 1.43(1) 1.75(1) 1.70(5) 0.99(2) 0.50(1) 



TABLE IL Critical exponents for directed sandpiles with bulk dissipation in d = 2. DR: 
Dhar and Ramaswamy's exact solution; DDS, deterministic directed model; ESDS and NESDS, 
stochastic directed models. Figures in parenthesis denote statistical uncertainties. 
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without spurious effects due to the finite lattice size. 

The moment analysis can be straightforwardly generalized to systems with bulk dissi- 
pation. In this case the role of the system size L as scaling parameter is played by the 
dissipation e. If the FSS hypothesis holds, the g-th moment for, say the size distribution, 
has an explicit dependence on the dissipation rate that reads: 

~ ^ ^-p.W. (10) 

The new moment Ps{q) = ^s{q + 1 ^ t~s) can be estimated by linear regression in a log- 
log plot of (s*^)^ as a function of e^^. Once this moment is computed, the exponent is 
given by = dps{q)/dq. The relation (s) = imposes Ps(l) = 1, and from here, once 
known A^, we compute Tg using the relation Ps(l) = A,(2 — r^). Analogous considerations 
allow to compute the exponents of the time distribution, A^ and r^. Finally, to check the 
exponents with the data collapse technique, one must plot the rescaled functions P(s)e~^''^'' 
as a function of s/e"^" and P(t)e~^*'^* as a function of t/e~^\ respectively. 
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In Table |T| we report the exponents computed in = 2 for the directed models DDS, 
ESDS, and NESDS with bulk dissipation. The corresponding moments Ps^q) and pt{<i) are 
shown in Figures]^, while Figs. |^ and ^ plot the data collapse for sizes and times, respectively. 

To conclude our analysis of directed sandpiles with bulk dissipation, we have proceeded 
to compute the correlation length of the models. In this case, the scaling of the correlation 
lengths with vanishing dissipation define the scaling exponents 

e||~e-^i, U-e-^'^- (11) 

and an affinity exponent Q = Using an analogous definition as in the case of boundary 

dissipation, we compute the exponents z/y = 1, z/^ = 1/2, and C, = 1/2, as shown in 
Fig. p. That is, the correlation length exponents are identical for both boundary and bulk 
dissipation. These results again imply an affinity exponent ^ = 1/2 in all the models studied 
so far. 

In view of these results, we have confirmation that the stochastic models belong to 
a different universality class than the deterministic directed sandpile. These results also 
point out in a very clear way that the critical behavior of models with boundary or bulk 
dissipation is identical. In fact, all critical exponents Ts,Tt,z, and ( are identical in both 
cases [0]. This further confirms the complete equivalence of both points of view with respect 
to sandpiles and shows that, at least in the directed case, the open boundary conditions 
usually implemented in simulations do not affect the scaling behavior in a peculiar way. 
Of course, the open boundary conditions breaks the translational invariance of the system, 
but in the thermodynamic limit this effect is negligible for the asymptotic critical behavior. 
Finally, these results validate theoretical approaches in which it is assumed a homogeneous 
dissipation that is much easier to treat analytically. 

As a last observation it is worth remarking that also in this case, a series of exponents 
such as ( and z/^ assume values independently of the universality class of the model under 
study. This sort of super-universality can be explained in terms of energy conservation as 



we shall see in Sec. VI 
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V. NUMERICAL SIMULATIONS OF ANISOTROPIC MODELS 

An important question to study in directed sandpile models is the effect on the scaling 
properties of any amount of diffusion along the preferred direction of transport x\\. One 
would expect that the broken symmetry introduced by the preferential direction should 
prevail on large scales, so that the dynamical scaling in directed and simply anisotropic 
sandpiles become indistinguishable in the thermodynamic limit. This fact hints towards the 
possibility of a unique universality class for both directed and anisotropic sandpiles. This 
universality class is determined uniquely by the lack of symmetry along the xy direction, 
and the presence or absence of stochastic elements in the definition of the models. 

In order to test this conjecture, we have performed numerical simulations of an 
anisotropic stochastic sandpile model, defined according to the following rules: on a hyper- 
cubic lattice of size L, we consider a model with threshold Zc = 2. When a site topples, it 
sends two grains of energy to two sites, randomly selected among the 2d + 1 nearest and 
next-nearest neighbors on the hyper-plane x\\ + 1, and the nearest neighbor on the hyper- 



plane X|| — 1, see Fig. [ly. The rules in this model are defined non-exclusive, in such a way 
that the same site can receive the two sand grain expelled by an active site. The model is 
clearly anisotropic, because the probability to transfer energy in the downwards direction is 
three times larger that in the upwards direction. It would thus correspond to a non-exclusive 
stochastic anisotropic sandpile (NESAS). We consider only the case of boundary dissipation, 
performing simulations for sizes ranging from L = 100 up to 6400, and averaging over 10^ 
nonzero avalanches. 

In Fig. |ll| we plot the correlation length and measured according to the rules 
given in Eqs. (^. We confirm the expectation that anisotropic models have the same scaling 
properties, as regards the scaling of the correlation lengths, as directed models with the same 
deterministic or stochastic ingredients. We have also measured the exponents Ts,Tt,D, and 
z for this model, using the moment analysis technique. The values found are = 1.43(1), 
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D = 1.75(1), Tt = 1.72(2), z = 0.98(2). These results, compared with Tables | and ||, show 
that this anisotropic models belongs to the same universality class of the ESDS and NESDS 
directed models, confirming the irrelevance of the diffusion along the preferred direction x» . 



VI. THE ROLE OF CONSERVATION IN SANDPILE MODELS 

We have seen in the preceding Sections that a subset of critical exponents characterizing 
the critical behavior of directed and anisotropic models have an interesting super-universal 
property; i.e. they are independent of the universality class of the models. In order to 
understand this feature we perform a theoretical analysis based on the conservation of energy, 
that is the basic symmetry in standard sandpile automata. We shall see in the following 
that the super- universal character of some critical exponents is dictated by simple energy 
conservation considerations. The use of this approach allows also to establish a relation 
between boundary and bulk dissipation models by introducing an effective dissipation that 
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depends on the system size. 

The avalanche dynamics in sandpile models is implicitly due to the imposed infinite time 
scale separation between driving and dissipation [p6| , p7| , p5| . In order to devise a theory 
that can take into account the symmetry introduced by the energy conservation, one must 
first regularize the rules of the models in such a way that a single time scale is ruling the 
dynamics. One way to do so is to introduce a nonzero driving rate, defined as the probability 
per unit time /i of a site to receive a grain of energy This driving rate plays the 

role of an external field and leads to the SOC behavior in the limit /i — *• 0+. On the other 
hand, given that the toppling rules are conserved, energy can leave the system only at the 
boundaries. Boundary dissipation is a natural choice in computer simulations. However, it 
introduces undesirable complications due to its singular character in a local theory. It is 
therefore convenient to use an homogeneous effective dissipation e, defined as the average 
energy lost in each toppling event. As observed in previous Sections, one can define models 
with periodic boundary conditions and built-in bulk dissipation. When constructing the 
local theory for models with open boundary conditions, the bulk dissipation e amounts to 
an effective parameter that is to be related to the system size L. 

With all these ingredients, we are ready to formulate conservation of energy as a con- 
tinuous equation. In sandpiles, we define the order parameter pa as the density of active 
sites (i.e., whose height z > z^). The only dynamics in the model is obviously due to the 
field pa{x,t), which is coupled to the local energy density, E{x,t) (i.e. the local density of 
sandgrains), which enhances or suppresses the generation of new active sites. A Langevin 
description for sandpile automata is possible by considering the dynamics of the local order- 
parameter field pa{x,t) in a coarse-grained picture, bearing in mind that the energy density 



E{x,t) is a conserved field. In Refs. |27,28], in analogy with absorbing-state phase transi- 
tions P^ , |57[ |, a pair of coupled dynamical equations for the fields pa{x,t) and E{x,t) were 
proposed. In the following we elucidate the consequences of energy conservation and we 
focus only on the latter equation. The interested reader can find the full set of equations 
in Ref. [^]. In the next subsections we shall consider separately directed and anisotropic 
models. 



A. Directed sandpiles 

We seek a continuous equation for the coarse-grained local density of energy, E{x,t). 
In the limit of zero driving and dissipation, energy is conserved. Therefore, the evolution 
equation fulfilled by the local field E is: 

^E^H = _v ■ - epaix, t) + h{x, t) + r]E{x, t). (12) 

The first term simply represents the diffusion of energy; the second term accounts for the 
dissipation that is associated with every toppling event; the third term represents the ex- 
ternal driving. Finally, the last term is a source of stochastic noise, that accounts for the 
randomness in the flow of energy. The noise term can be generated by the toppling rules in 
a stochastic model, or by the initial conditions plus the random driving in a deterministic 
model. We will require the noise to have zero average: 
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(r/s(f,t)) = 0. 



(13) 



The noise correlator {riE{x,t)r]E{x' ,t') is of fundamental importance for the determination 
of universality classes and the critical behavior of the order parameter. However, for our 
present purposes we do not need precise knowledge of its analytical form (for a detailed 



discussion, see Refs. ||27| , |28[| ) . 

The current can be constructed by appealing to the symmetries of the model. The 
transport of energy is due to topplings. These are isotropic along the transversal direction 
x±, therefore the current along this direction will be proportional to the gradient of the 
density of active sites. In the preferred direction, on the other hand, all the energy is 
transferred downwards; therefore, the current in this direction must be proportional to the 
density of active sites. The final form of the current is then 

Js(f,t) = -D±V±pa{x,t) + 2Xpa{x,t)e\\. (14) 

Plugging this expression into the equation for the energy, we have the final result: 

dE{x, t) 



dt 



D^V\pa{x, t) - 2\d\\Pa{x, t) - epa{x, t) + h{x, t) + 7]e{x, t) , (15) 



where the symbol d\\ stands for the partial derivative d/dx\\. This is the general conservation 
equation for any directed sandpile model. It is worth remarking at this point that the energy 
field is a static field, in the sense that energy diffuses only if active sites are present in the 
system. This is intuitively understood in sandpile models, where energy (sand) grains diffuse 
only from toppling sites. 

To analyze the consequences of Eq. (0), it proves useful to define the susceptibility 

where the symbol ()^ denotes an average over the noise distribution. By definition, the 
susceptibility measures the average increase in the number of active sites due to an impulsive 
perturbation, that is, to the addition of a single energy grain. Since we measure the size of 
the avalanches by the total number of topplings, the average avalanche size is given by 

(s) = / d'^xdtxix,t). (17) 



Taking the functional derivative of Eq. (|l^) and averaging over time and noise, we obtain, 
in the limit t — oo, in which the sandpile is in a stationary state with constant average 
energy, the following equation for the static susceptibility: 

D,_Vlx{S) - 2A9||x(x) - ex{x) = -6^''\x). (18) 

This equation can be easily solved in Fourier space. Defining the transformation 

Xixii^x^) = J d''-'kdqx{q,k)e''-'^e''^^^l (19) 
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and substituting into Eq. (0), we obtain the solution 



Xil,k)= ^ , ^ , (20) 



D±k^ + 2iXq + e 
which yields the susceptibility in real space 

X(H ■ = / ^^-fe"-^^ ^^^J^T^^, (21) 

The integral in q is easily performed by residues. The remaining d—1 integrals in k become 
then decoupled Gaussian integrals that yield the result, setting = 1: 

-I / X \ {d-l)/2 

Xix^uX^) = ^^] 3.(l-'^)/2e-X|,./2Ag-Axi/2.|,_ (22) 

2A yzTT J " 
Eq. can be conveniently rewritten into the scaling form 

x(xi|,f.)=xf-^)/^r(|,|^), (23) 

where F is a cut-off function that decreases exponentially in both its arguments. Comparing 
this last expression with Eq. (P^), we can identify the parallel and transversal correlation 
lengths: 

~e-\ (24) 



In more general terms, if we define the exponents z/y and by Eqs. (|TT]), then we have for 
directed sandpiles i/y = 1 and = 1/2. From these last expressions, we can read off a first 
exact result for directed sandpiles: the avalanches produced in those models are elongated, 
with characteristic length in the parallel and transversal directions related by an affinity 
exponent ( = 1/2. It is very important to stress that these results are independent of the 
particular model consider and of the dimensionality d of the system, dictated only by the 
energy balance in the stationary state. 



We can use the result ([2^ ) to relate the effective bulk dissipation with the system size 
in a model with open boundary conditions. To sustain a steady state with constant average 
energy, avalanches must reach the bottom boundary in order to be able to dissipate. This 
means that the characteristic length of the avalanches in the parallel direction must be 
proportional to the system size ~ L. We have therefore that in boundary dissipation 
models we can define an effective dissipation rate e that is related with the system size by: 

e ~ L"\ (25) 

From this relations we easily find that Ag = D and = z. These identities are recovered 
in numerical simulations (see Tables | and |I|). Finally, from Eq. (pO|), we can recover the 
well-known result linking the system size and the average avalanche size, (s) = x{q = 0,k = 



0) ~ e~ ~ L [10,15,23|. That is, for any directed model, the average avalanche size must 
scale as the inverse of the dissipation, in the case of model with bulk dissipation, or as 
the size of the system along the preferred direction of transport for sandpiles defined with 
boundary dissipation. 
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B. Anisotropic sandpiles 



Having completed the analysis of directed sandpiles, we turn our attention to the more 
complex case of anisotropic sandpiles. In this kind of model, the transport of energy is not 
strictly directed in the parallel direction, but is simply stronger in the direction than in 
the opposite direction —x\\ . The presence of backwards flow allows the possibihty of diffusion 
in the preferred direction, and thus the equation for the conservation of energy becomes in 
this case 



dE{x,t) 
di 



D^Vlpa{x, t) + D\^dlpa{x, t) - 2Xd\\Pa{x, t) - epa{x, t) + h{x, t) + //^(f , t). (26) 



Even though it is straightforward to obtain Eq. pq) by symmetry arguments, it is instructive 
to derive it also starting from the anisotropy of the microscopic rules. To this end, let us 
consider a one dimensional model in which active sites, when toppling, send a fraction of 
energy p to the left neighbor, and a fraction (1 — p) to the right neighbor, with < p < 1. 
Let us coarse-grain the line into cells of a certain (small) size. The variation of energy of 
the i-ih cell after a parallel updating will be given by 



-f^ + (1 



e)pp^+^) + {l-e){l-p)p^: 



(i-i) 



(27) 



where p^'-* is the density of active sites in the i-th cell. The proportionality stems from the 
fact that the actual contribution comes only from the boundary sites in the cell. This last 
equation can be rewritten 



-epf - (1 - e)p 



ft'' 



(28) 



For p = or p = 1 (for a strictly directed model), only one term remains and, after passing 
to a continuous formulation, we recover Eq. (|15|), restricted to the preferred direction. For 
p 7^ 0, 1, and after performing some algebraic manipulations, we obtain 



AE(') - -epf + (1 
- (1 - e)p 



p« 



+ (1 



2P« 
e)(l-p) 

ii) 



6pW + (1 _ e)p9^pW + (1 _ + (1 - e)(l - p){-d^p^: 



epf + (1 - e)pdlpt^ - (1 - e)(l - 2p)d^pt 



(i) 



where in the third line we have passed to the continuum limit, defining the discretized first 



derivative dxP^a 



P«- 



-P'a 



1*+^) and second derivative d^p^^^ 



-Pt'^- 



2p«. Forp = 1/2 



we recover, as expected, an isotropic equation. For all other values of p, the equation for 
the energy will be effectively given by (|26|) , with the phenomenological parameter A given 
at the microscopic level by the expression (1 — e) (1/2 — p). 

From Eq. (^), we can obtain the corresponding equation for the susceptibility. The 
solution in Fourier space is readily found to be 



1 



(29) 



17 



Upon integration over k and q, one obtains the expression in real space 

X(xiux^) = ——: / d'^'-^fce*"^ / dq ^ . (30) 

This last integral can be performed analytically in = 1 and 2 (see Appendix). For d > 2, 
even though we do not have a closed expression, we can obtain the leading scaling behavior. 
To simplify the calculations, we set, without lack of generality, D± = D\\ = 1. The integra- 
tion in q is done by the method of the residues (see Appendix). The integration of the k 
angular part |Q yields 

X{x\\,x±) = -[ — ) xl" dzz^^^J^i-fXxz) , . (31) 



2 \2nJ ^ Jo -^-K'-^-J (1 + ^2)1/2 

Here, Ju{z) is the first kind Bessel function of order z/, and we have defined the constants 
1/ = {d — 3)/2 and 7 = (A^ + e)^/^. We are interested in the behavior of this integral for 
large distances, that is, in the limit x\\ ^ x_l ^ 1. In this limit, the weight of the integral is 
given by the region of small z, since the exponential suppresses large values. We can then 
approximate the integral in the interval < 2; < 1 and perform a Taylor expansion of the 
square root in the exponential and the denominator. In the denominator, we readily have 
(1 + z'^y^'^ ~ 1. The term in the exponential, however, contains a constant term, and must 
be therefore expanded up to second order: 

- X|i (7VI + z^ -X)- (7[1 + - A) = (7 - A) - X|,7zV2. (32) 

In the limit e — * 0, we have 7 ~ A, and the constant 7 — A can be expanded to give: 

.-A = (A= + .)V^-A.a(i + ^)-A=|^. (33) 

Substituting these approximations into Eq. (|3l|), we are led to the expression 



x(x||,x^) ^ T^TTr^TT^l'"'^^"'"^^'' Hdy y'^'My)e'^-'\ (34) 



2A (27r) 



where we have performed the change of variables y = 7x^2; and extended again the upper 
limit of the integral to infinity (which is allowed given its exponential convergence). The 
integral in Eq. ( |3lD yields p8[] : 

which as usual we can write in the scaling form, 

/ Tit nr 1 \ 

(36) 



From here, the correlation exponents read z/y = 1 and z/^ = 1/2, as in the directed case. 
This implies again an affinity exponent ( = 1/2. 
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FIG. 12. Color plots of the local density of topplings in two avalanches of size 50 000 for the 
(a) DDS and (b) ESDS models. White stands for a sigle toppling per site; black represents the 
maximum number of topplings. 



The conclusion of the lengthy calculations developed in this section is that the presence 
of any amount of diffusion along the preferred direction of a directed sandpile model is 
completely irrelevant. As soon as there is anisotropy in a model (in our mathematical for- 
mulation, when A 7^ 0, however small), it takes over and places the model in the universality 
class of completely directed sandpiles. In particular, we recover the result (s) ~ L for any 
anisotropic sandpile, in agreement with the numerical results in Ref. [l^ and the analytic 
results of Ref. pS]. We remark, however, that our results do not rely in a particular model 



like those of |23|, but only on symmetry arguments, and are therefore of a broader generality. 



VII. CONCLUSIONS 

In this paper we have presented a detailed numerical analysis of deterministic and 
stochastic directed sandpile models. We find definitive evidence for the existence of a new 
universality class, embracing directed sandpile models with stochastic rules. The origin of 
the different critical behavior can be traced back to the presence of multiple topplings in 



the latter case. An example of this feature is provided in Fig. |T2|, where we plot the local 
density of topplings in two avalanches corresponding to the DDS and ESDS models. From 
this figure it becomes evident that the stochastic dynamics induces multiple toppling events, 
which are forbidden in the deterministic models. This feature has been fruitfully exploited 
in Ref. EM to obtain an analytical solution of the stochastic model. 
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We have also studied the case of directed sandpiles with bulk dissipation. In this case, 
our results prove that the critical behavior is unchanged. This points out that the boundary 
dissipation does not play any particular role in the development of the critical behavior in 
directed sandpiles. 

Finally, numerical results indicate that some critical exponents show a super-universal 
nature, assuming the same values independently of the universality class. We provide an 
analytical explanation of this feature by means of a continuous phenomenological equation 
that takes into account the energy balance condition imposed by the dynamical rules in 
sandpile models. 
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APPENDIX: 

In this Appendix, we work out the exact value of the integral (pO[), giving the static 
susceptibility for anisotropic sandpiles, in ci = 1 and 2. We take again, without lack of 
generality, = D\\ = 1. 



Static susceptibility in d = 1 

In (i = 1, Eq. (pOD reduces to the simple form 

X{x\\) = — / dq— ^ . (Al) 

This integral can be done by the method of the residues. The zeroes of the denominator are 

g± = ^[±(A2 + e)i/2-A]. (A2) 
Closing the contour of the integral arround the upper complex plane, we obtain 



1 -X|i{VA2+e-A) 

X(X||) = - , (A3) 



We are interested in the limit e ^ 0. The denominator gives a/A^ + e ^ A, while the 
argument in the exponential yields 



VA2+7-A = A 

From here, we write 



e \ V2 



2A- <*^) 



(A5) 



which allows to determine the correlation length ^|| = e ^• 
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Static susceptibility in d = 2 

In the case = 2, we have 



1 poo roo p^Q^W 



Integrating by residues the q variable, similarly as in the previous subsection, we obtain 

y(xii,x_l) = — e^""'" / dk cos(kx±) , . (A7) 

The last integral is duable giving a modified Bessel function of order 0: 

= ^e^^"^o [(A^ + 6)^/^(4 + . (A8) 

In the limit ^ a;_L ^1, the modified Bessel function Kq can be replaced by its asymptotic 
form Kq(z) ^ e"^/(27r2;)^/2 ||33]. Then, we have 



^(^11'^^) ^ (2703/2 (A2 + e)i/4 + 4)1/4 (^^H " f^' + e] + 4] V^) . (A9) 

In the double limit x\\ ^ x± ^ 1 and e ^ 0, the argument in the last exponential can be 
expanded 

A., - [A' + .ri4 . 41"= - A^. - (a ^) + ||) - (M") 

where we have only kept terms linear in e and x\/x\\. 
The static susceptibility can be finally written 

x(-lh-x) - (^A^-F'''-"^"^'"-""^^'""' (All) 
from which we immediately read the correlation lengths .^n = and = e^^/^ 
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